# Replication File for "Preaching Politics: How Politicization Undermines Religious Authority in the Middle East"
# Updated April 2022
# Contact Scott Williamson with questions: scott.williamson@unibocconi.it

### Install Packages:
rm(list = ls())
library(tidyverse)
library(lmtest)
library(sandwich)
library(ggplot2)
library(ggthemes)
library(ggpubr)
library(stargazer)
library(cregg)
library(dotwhisker)
library(dplyr)

### Set Working Directory
# setwd()

source("functions/vcovCluster.R")
source("functions/effect_plotter.R")

### Load Data
conjoint = readRDS("religion_conjoint_reshaped.rds")

head(conjoint)

### Cleaning Attributes

conjoint = conjoint %>% 
  mutate(trust = recode(trust, "Not trust at all" = "0", "A little" = "1", "Somewhat" = "2", "A lot" = "3",
                        "Completely trust" = "4", .default = NA_character_),
         trust = as.numeric(trust),
         cjt_nationality = recode(cjt_nationality, "$q_country" = "Home"),
         cjt_age = factor(cjt_age, levels = c("28", "45", "70")),
         cjt_education = factor(cjt_education, levels = c("Self-taught in religion", "BA in Engineering", "BA in Islamic Fiqh",
                                                          "Doctorate in Engineering", "Doctorate in Islamic Fiqh")),
         cjt_nationality = factor(cjt_nationality, levels = c("Home", "Iraq", "Saudi Arabia")),
         cjt_hafiz = factor(cjt_hafiz, levels = c("No", "Yes")),
         cjt_location = factor(cjt_location, levels = c("Internet", "Government Television", "Satellite Television",
                                                        "Mosque registered with government", "Mosque not registered with government")),
         cjt_followers = factor(cjt_followers, levels = c("5,000", "50,000", "500,000")),
         cjt_ideology = factor(cjt_ideology, levels = c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                        "Sufi", "Shiite")),
         cjt_usa = factor(cjt_usa, levels = c("No Opinion", "Participates in meetings with US embassy",
                                              "Criticizes United States", "Supports armed resistance against United States")),
         cjt_politics = factor(cjt_politics, levels = c("None", "Supports Government", "Criticizes Government")))


age = unique(conjoint$cjt_age)
education = unique(conjoint$cjt_education)
nationality = unique(conjoint$cjt_nationality)
hafiz = unique(conjoint$cjt_hafiz)
location = unique(conjoint$cjt_hafiz)
followers = unique(conjoint$cjt_followers)
ideology = unique(conjoint$cjt_ideology)
usa = unique(conjoint$cjt_usa)
politics = unique(conjoint$cjt_politics)


refcats = data_frame(variable = c("28", "Self-taught in religion", "Home", "No", "Internet", "5,000",
                                  "Sunni Orthodox", "No Opinion", "None"),
                     coef = 0, se = 0, attribute = 0)

attributes = data_frame(variable = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                     "Followers", "Ideology", "USA", "Politics"),
                        coef = NA_real_,
                        se = NA_real_, attribute = 1)
  
# Subset data to Sunni Muslims
conjoint = conjoint %>% filter(sect == "Sunni"|sect=="Salafi")

### Primary Conjoint Results: Figure 1 from Article

df = conjoint

# Regression Model
basemod1 = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
mod1<-summary(basemod1)
basemod1 = coeftest(basemod1, vcov = vcovCluster(basemod1, factor(df$id)))
basemod1

# Robustness Check with Fixed Effects for Respondent Country
fe_mod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                         cjt_followers + cjt_ideology + cjt_usa + cjt_politics + factor(residence_country), data = df)
fe_mod = coeftest(fe_mod, vcov = vcovCluster(fe_mod, factor(df$id)))
fe_mod
stargazer(fe_mod)

# Plot
p1 = effect_plotter(model.output =  basemod1, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Forced Choice Conjoint Results")
p1
ggsave(filename="figs_resubmission/conjoint_main.jpg",p1,scale = 1, width = 7, height = 8, units = "in", dpi = 500)

# Marginal Means for Primary Results

basemod1_mm<- binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                cjt_followers + cjt_ideology + cjt_usa + cjt_politics

p1_mm<-plot(mm(df, basemod1_mm, id = ~ id), vline = 0.5)
p1_mm
ggsave(filename="figs_resubmission/conjoint_main_mm.jpg",p1_mm,scale = 1, width = 7, height = 8, units = "in", dpi = 500)



### Subgroup Effects for Islamists and Non-Islamists 

# Subsetting to Islamists
df = conjoint %>% filter(islamist==1)

# Primary Results Among Islamists (coefficients reported in Figure 2 of Article)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))
basemod

plot_islamist = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Islamists")


# Subsetting to Non-Islamists

df = conjoint %>% filter(islamist==0)

# Primary Results Among Non-Islamists (coefficients reported in Figure 2 of Article)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))
basemod

plot_not_islamist = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = F, title = "Non-Islamists")
#ggarrange(p1, p2, widths = c(3.55, 1.45))%>% 
islamist_plots<-ggarrange(plot_islamist,plot_not_islamist,ncol=2,widths = c(3.05, 1.6))
ggsave("figs_resubmission/islamist__plots.jpg",islamist_plots, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

# Figure 2 from Main Article

term<-c("Sunni Orthodox","Muslim Brotherhood","Salafi","Sufi","Shiite","Sunni Orthodox","Muslim Brotherhood","Salafi","Sufi","Shiite")
estimate<-c(0,-0.094,-0.109,-0.203,-0.396,0,-0.161,-0.124,-0.140,-0.343)
std.error<-c(0,0.014,0.014,0.014,0.013,0,0.009,0.009,0.009,0.009)
model<-c("Islamist","Islamist","Islamist","Islamist","Islamist","Non-Islamist","Non-Islamist","Non-Islamist","Non-Islamist","Non-Islamist")

islamist_df<-cbind.data.frame(term,estimate,std.error,model)
islamist_df

islamist_df$term<-factor(islamist_df$term,levels=c("Sunni Orthodox","Muslim Brotherhood","Salafi","Sufi","Shiite"))

islamist_subgroup_plot<-dwplot(islamist_df, dot_args = list(size = 2.5),whisker_args = list(size = 1)) + 
  theme_bw() + #geom_point(size = 3, shape = 16) +
  geom_vline(xintercept = 0, colour = "black", linetype = 2) +
  ggtitle("") + 
  scale_x_continuous("AMCE",limits = c(-.45, .05)) + theme(panel.spacing = unit(1, "lines")) +
  scale_y_discrete("Ideology Levels\n",labels=c("Shiite","Sufi","Salafi","Muslim \nBrotherhood","Sunni \nOrthodox")) +
  theme(axis.title=element_text(size=10),
        legend.text=element_text(size=10),legend.position="bottom",
        axis.text.x=element_text(size=10,color="black"),
        axis.text.y=element_text(size=10,color="black")) + scale_color_manual(values=c("darkgrey", "black")) +
  labs(color='Respondent Ideology')
islamist_subgroup_plot
ggsave("figs_resubmission/islamist_subgroup_plot.jpg",islamist_subgroup_plot, scale = 1, width = 7, height = 5, units = "in", dpi = 500)

# Marginal Means for Islamist Plot

conjoint$islamist_factor<-factor(conjoint$islamist)
islamist_mm<- cj(conjoint,binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
  cjt_followers + cjt_ideology + cjt_usa + cjt_politics,estimate="mm",by=~islamist_factor)

islamist_plots_mm<-plot(islamist_mm,group="islamist_factor", vline = 0.5)
islamist_plots_mm
ggsave("figs_resubmission/islamist__plots_mm.jpg",islamist_plots_mm, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

# Regression Model for Islamist Interaction

islamist_mod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                    cjt_followers + cjt_ideology*islamist + cjt_usa + cjt_politics, data = conjoint)
islamist_mod = coeftest(islamist_mod, vcov = vcovCluster(islamist_mod, factor(conjoint$id)))
islamist_mod
stargazer(islamist_mod)
  


### Subgroup Effects for Government Opponents and Supporters

# Subsetting to Opponents

# Primary Results Among Opponents (coefficients reported in Figure 3 of Article)
df = conjoint %>% filter(critic==1)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))
basemod

plot_critic = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Government Critics")

# Subsetting to Supporters
df = conjoint %>% filter(critic == 0)

# Primary Results Among Supporters (coefficients reported in Figure 3 of Article)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))
basemod

plot_support = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = F, title = "Government Supporters ")
#ggarrange(p1, p2, widths = c(3.55, 1.45))%>% 
# ggsave(filename = "figs_scott/conjoint_govt.png")
gov_plots<-ggarrange(plot_critic,plot_support,ncol=2,widths = c(3.05, 1.6))
ggsave("figs_resubmission/gov__plots.jpg",gov_plots, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

# Figure 3 from Main Article

term<-c("None","Supports Government","Criticizes Government","None","Supports Government","Criticizes Government")
estimate<-c(0,-0.045,-0.008,0,-0.003,-0.062)
std.error<-c(0,0.009,0.008,0,0.010,0.010)
model<-c("Opponent","Opponent","Opponent","Supporter","Supporter","Supporter")

politics_df<-cbind.data.frame(term,estimate,std.error,model)
politics_df

politics_df$term<-factor(politics_df$term,levels=c("None","Supports Government","Criticizes Government"))

politics_subgroup_plot<-dwplot(politics_df, dot_args = list(size = 2.5),whisker_args = list(size = 1)) + 
  theme_bw() + #geom_point(size = 3, shape = 16) +
  geom_vline(xintercept = 0, colour = "black", linetype = 2) +
  ggtitle("") + 
  scale_x_continuous("AMCE",limits = c(-.1, .05)) + theme(panel.spacing = unit(1, "lines")) +
  scale_y_discrete("Political Opinion Levels\n",labels=c("Criticizes \nGovernment","Supports \nGovernment","None")) +
  theme(axis.title=element_text(size=10),
        legend.text=element_text(size=10),legend.position="bottom",
        axis.text.x=element_text(size=10,color="black"),
        axis.text.y=element_text(size=10,color="black")) + scale_color_manual(values=c("darkgrey", "black")) +
  labs(color='Respondent Support for Government')
politics_subgroup_plot
ggsave("figs_resubmission/politics_subgroup_plot.jpg",politics_subgroup_plot, scale = 1, width = 7, height = 5, units = "in", dpi = 500)

# Marginal Means for Gov Plot

conjoint$critic_factor<-factor(conjoint$critic)
critic_mm<- cj(conjoint,binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                   cjt_followers + cjt_ideology + cjt_usa + cjt_politics,estimate="mm",by=~critic_factor)

critic_plots_mm<-plot(critic_mm,group="critic_factor", vline = 0.5)
critic_plots_mm
ggsave("figs_resubmission/critic__plots_mm.jpg",critic_plots_mm, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

# Regression Model for Views of Government Interaction

df = conjoint %>% filter(critic==1|critic==0)
critic_mod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                  cjt_followers + cjt_ideology + cjt_usa + cjt_politics*critic, data = df)
critic_mod = coeftest(critic_mod, vcov = vcovCluster(critic_mod, factor(df$id)))
critic_mod
stargazer(critic_mod)


### Subgroup Effects for Views of Terrorism

# Subsetting to Terrorism Sometimes Justified

df = conjoint %>% filter(terror_never==0)

# Primary Results Among Terrorism Justified (coefficients reported in Figure 4 of Article)

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))
basemod

plot_justified = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Terrorism Justified")

#Subsetting to Terrorism Never justified

df = conjoint %>% filter(terror_never == 1)

# Primary Results Among Terrorism Never Justified (coefficients reported in Figure 4 of Article)

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))
basemod

plot_never = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = F, title = "Terrorism Never Justified")
#ggarrange(p1, p2, widths = c(3.55, 1.45))%>% 
#  ggsave(filename = "figs_scott/conjoint_terrorism.png")
terror_plots<-ggarrange(plot_justified,plot_never,ncol=2,widths = c(3.05, 1.6))
ggsave("figs_resubmission/terror__plots.jpg",terror_plots, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

# Figure 4 from Main Article

term<-c("No Opinion","Participates in meetings with US embassy","Criticizes United States","Supports armed resistance against United States","No Opinion","Participates in meetings with US embassy","Criticizes United States","Supports armed resistance against United States")
estimate<-c(0,-0.020,-0.004,-0.038,0,-0.028,-0.014,-0.097)
std.error<-c(0,0.011,0.011,0.011,0,0.008,0.008,0.009)
model<-c("Sometimes Justified","Sometimes Justified","Sometimes Justified","Sometimes Justified","Never Justified","Never Justified","Never Justified","Never Justified")

terror_df<-cbind.data.frame(term,estimate,std.error,model)
terror_df

terror_df$term<-factor(terror_df$term,levels=c("No Opinion","Participates in meetings with US embassy","Criticizes United States","Supports armed resistance against United States"))

terror_subgroup_plot<-dwplot(terror_df, dot_args = list(size = 2.5),whisker_args = list(size = 1)) + 
  theme_bw() + #geom_point(size = 3, shape = 16) +
  geom_vline(xintercept = 0, colour = "black", linetype = 2) +
  ggtitle("") + 
  scale_x_continuous("AMCE",limits = c(-.15, .05)) + theme(panel.spacing = unit(1, "lines")) +
  scale_y_discrete("USA Levels\n",labels=c("Supports armed \nresistance \nagainst USA","Criticizes USA","Participates in \nmeetings with \nUS embassy","No Opinion")) +
  theme(axis.title=element_text(size=10),
        legend.text=element_text(size=10),legend.position="bottom",
        axis.text.x=element_text(size=10,color="black"),
        axis.text.y=element_text(size=10,color="black")) + scale_color_manual(values=c("darkgrey", "black")) +
  labs(color='Respondent Views of Terrorism')
terror_subgroup_plot
ggsave("figs_resubmission/terror_subgroup_plot.jpg",terror_subgroup_plot, scale = 1, width = 7, height = 5, units = "in", dpi = 500)

# Marginal Means for Terrorism Plot

conjoint$terror_never_factor<-factor(conjoint$terror_never)
terror_mm<- cj(conjoint,binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                   cjt_followers + cjt_ideology + cjt_usa + cjt_politics,estimate="mm",by=~terror_never_factor)

terror_plots_mm<-plot(terror_mm,group="terror_never_factor", vline = 0.5)
terror_plots_mm
ggsave("figs_resubmission/terror__plots_mm.jpg",terror_plots_mm, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

# Regression Model for Views of Terrorism

terror_mod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                  cjt_followers + cjt_ideology + cjt_usa*terror_never + cjt_politics, data = conjoint)
terror_mod = coeftest(terror_mod, vcov = vcovCluster(terror_mod, factor(conjoint$id)))
terror_mod
stargazer(terror_mod)



#####################
# Robustness Checkes#
#####################

### Results By Education Level

# Subsetting to University Graduates

df = conjoint %>% filter(college==1)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))


p1 = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "College Degree")

# Subsetting to Non-University Graduates

df = conjoint %>% filter(college == 0)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))


p2 = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = F, title = "No College Degree")
edu_plots<-ggarrange(p1, p2, widths = c(3.30, 1.7))
ggsave("figs_resubmission/edu_subgroup.jpg",edu_plots, scale = 1, width = 7, height = 8, units = "in", dpi = 500)


### Results By Gender

# Subsetting to Men

df = conjoint %>% filter(men==1)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))


p1 = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Men")

# Subsetting to Women 

df = conjoint %>% filter(men == 0)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))


p2 = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = F, title = "Women")
gender_plots<-ggarrange(p1, p2, widths = c(3.30, 1.7)) 
ggsave("figs_resubmission/gender_subgroup.jpg",gender_plots, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

### Results By Religiosity

# Subsetting to Religious Respondents

df = conjoint %>% filter(rlgn_pca > 0)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))


p1 = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "More Religious")

# Subsetting to Non-Religious Respondents

df = conjoint %>% filter(rlgn_pca < 0)
basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))


p2 = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = F, title = "Less religious")
relig_plots<-ggarrange(p1, p2, widths = c(3.30, 1.7)) 
ggsave("figs_resubmission/relig_subgroup.jpg",relig_plots, scale = 1, width = 7, height = 8, units = "in", dpi = 500)



### Results with Scale Outcome

df = conjoint

basemod2 = lm(trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
mod2<-summary(basemod2)
basemod2 = coeftest(basemod2, vcov = vcovCluster(basemod2, factor(df$id)))
basemod2

scalar_plot = effect_plotter(model.output =  basemod2, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -.75, x.upper = 0.5, labs = T, title = "Scalar Trust in Authority Outcome")
scalar_plot
ggsave("figs_resubmission/scalar_plot.jpg",scalar_plot, scale = 1, width = 7, height = 8, units = "in", dpi = 500)


### Results with Second Binary Outcome (Prefer to Hear Sermon)

df = conjoint

basemod3 = lm(binary_prefer ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
mod3<-summary(basemod3)
basemod3 = coeftest(basemod3, vcov = vcovCluster(basemod3, factor(df$id)))
basemod3

sermon_plot = effect_plotter(model.output =  basemod3, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Forced Choice Sermon Outcome")
sermon_plot
ggsave("figs_resubmission/sermon_plot.jpg",sermon_plot, scale = 1, width = 7, height = 8, units = "in", dpi = 500)


### Regression Table for all Three Outcomes (Reported in SI)

stargazer(basemod1, basemod2, basemod3,
          title = "Conjoint Regression Models",
          dep.var.labels   = c("Authority Outcomes"),
          covariate.labels = c("45", "70","Engineering BA","Islamic Fiqh BA","Engineering Doctorate",
                               "Islamic Fiqh Doctorate", "Hafiz", "Iraqi", "Saudi", "Government TV", "Satelitte TV",
                               "Mosque registered w/ Gov", "Mosque not registered w/ Gov", "50,000", "500,000", "Muslim Brotherhood",
                               "Salafi", "Sufi", "Shiite", "Meets with US Embassy", "Criticizes US", "Supports Armed Resistance Against US",
                               "Supports Gov", "Criticizes Gov"),
          font.size = "scriptsize",
          keep.stat = c("n","adj.rsq"),
          model.names = TRUE,
          multicolumn = TRUE,
          notes = "SEs clustered by survey respondent")


### Results from First Conjoint Table Only

df = conjoint[(conjoint$variable==1|conjoint$variable==2),]

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
                cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = df)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(df$id)))
basemod

first_conjoint_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                  "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Forced Choice Outcome")
first_conjoint_plot
ggsave("figs_resubmission/first_conjoint.jpg",first_conjoint_plot, scale = 1, width = 7, height = 8, units = "in", dpi = 500)

### Results by by Country

egypt_conjoint<-conjoint[conjoint$residence_country=="Egypt",]
jordan_conjoint<-conjoint[conjoint$residence_country=="Jordan",]
gulf_conjoint<-conjoint[conjoint$residence_country=="Bahrain"|conjoint$residence_country=="Kuwait"|conjoint$residence_country=="Qatar",]
morocco_conjoint<-conjoint[conjoint$residence_country=="Morocco",]
saudi_conjoint<-conjoint[conjoint$residence_country=="Saudi Arabia",]
tunisia_conjoint<-conjoint[conjoint$residence_country=="Tunisia",]
lebanon_conjoint<-conjoint[conjoint$residence_country=="Lebanon",]
turkey_conjoint<-conjoint[conjoint$residence_country=="Turkey",]
uae_conjoint<-conjoint[conjoint$residence_country=="United Arab Emirates",]

# Egypt

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = egypt_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(egypt_conjoint$id)))

egypt_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Egypt Results")
egypt_plot

# Jordan

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = jordan_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(jordan_conjoint$id)))

jordan_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Jordan Results")
jordan_plot

# Gulf Countries

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = gulf_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(gulf_conjoint$id)))

gulf_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.6, x.upper = 0.2, labs = T, title = "Gulf Results")
gulf_plot

# Morocco

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = morocco_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(morocco_conjoint$id)))

morocco_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Morocco Results")
morocco_plot

# Saudi Arabia

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = saudi_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(saudi_conjoint$id)))

saudi_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Saudi Results")
saudi_plot

# Tunisia

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = tunisia_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(tunisia_conjoint$id)))

tunisia_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Tunisia Results")
tunisia_plot

# Lebanon

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = lebanon_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(lebanon_conjoint$id)))

lebanon_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Lebanon Results")
lebanon_plot

# Turkey

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = turkey_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(turkey_conjoint$id)))

turkey_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "Turkey Results")
turkey_plot

# UAE

basemod = lm(binary_trust ~  cjt_age + cjt_education + cjt_hafiz + cjt_nationality + cjt_location + 
               cjt_followers + cjt_ideology + cjt_usa + cjt_politics, data = uae_conjoint)
basemod = coeftest(basemod, vcov = vcovCluster(basemod, factor(uae_conjoint$id)))

uae_plot = effect_plotter(model.output =  basemod, names.variables = c("Age", "Education", "Hafiz", "Nationality", "Location",
                                                                 "Followers", "Ideology", "USA", "Politics"), 
                    names.levels = list(c("28", "45", "70"), c("Self-taught in religion", "BA in Engineering", 
                                                               "BA in Islamic Fiqh","Doctorate in Engineering",
                                                               "Doctorate in Islamic Fiqh"),c("No", "Yes"), 
                                        c("Home", "Iraq", "Saudi Arabia"), 
                                        c("Internet", "Government Television", "Satellite Television",
                                          "Mosque registered w/ gov", "Mosque not registered w/ gov"),
                                        c("5,000", "50,000", "500,000"), c("Sunni Orthodox", "Muslim Brotherhood", "Salafi",
                                                                           "Sufi", "Shiite"),
                                        c("No opinion", "Meets with US embassy",
                                          "Criticizes US", "Supports armed resistance against US"),
                                        c("None", "Supports government", "Criticizes government")), 
                    effect.label = "Change in Probability", 
                    x.lower = -0.4, x.upper = 0.2, labs = T, title = "UAE Results")
uae_plot

# Combining Country Plots

country_plots1<-ggarrange(egypt_plot,jordan_plot,gulf_plot,morocco_plot,ncol=2,nrow=2)
country_plots1
ggsave("figs_resubmission/country__plots1.jpg",country_plots1, scale = 1, width = 10, height = 12, units = "in", dpi = 500)

country_plots2<-ggarrange(saudi_plot,tunisia_plot,lebanon_plot,turkey_plot,uae_plot,ncol=3,nrow=2)
country_plots2
ggsave("figs_resubmission/country__plots2.jpg",country_plots2, scale = 1, width = 14, height = 14, units = "in", dpi = 500)

